Overscreening in ID lattice Coulomb gas model of ionic liquids 
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Overscreening in the charge distribution of ionic liquids at electrified interfaces is shown to pro- 
ceed from purely electrostatic and steric interactions in an exactly soluble one dimensional lattice 
Coulomb gas model. Being not a mean- field effect, our results suggest that even in higher dimen- 
sional systems the overscreening could be accounted for by a more accurate treatment of the basic 
lattice Coulomb gas model, that goes beyond the mean field level of approximation, without any 
additional interactions. 



Room temperature ionic liquids (RTILs) are Coulomb 
fluids with large, asymmetric ions, and are used as elec- 
trolytes in fuel and solar cells, batteries and super capac- 
itors, to name but a few important applications [1]. As 
pointed out in several seminal contributions by Korny- 
shev [2] , the size of the ionic species leads in general to 
crowding and lattice saturation, thus engendering a fun- 
damentally different behavior of ionic liquids at charged 
interfaces as compared to aqueous electrolytes. In some 
particular cases several aspects of this behavior can be 
captured on a mean- field level by the lattice Coulomb gas 
(LCG) Kornyshev model 0,[3|. 

Steric effects stipulate that the capacitance of ionic liq- 
uids decays at large applied voltages while at the point of 
zero charge (PZC) it can exhibit a maximum as well as 
a minimum (and is thus a nonmonotonic function of ap- 
plied voltage) depending on the lattice packing fraction 
y , whereas for dilute electrolytes the capacitance at PZC 
is always a^minimum [2[. X-ray reflectivity [5], SFA [6| 
and AFM [7] studies at charged interfaces reveal an alter- 
nating charge distribution starting with an overscreening 
cat ionic layer, at the negatively charged substrate, which 
decays roughly exponentially into the bulk liquid with a 
periodicity comparable with the size of ionic species [8] . 
This observed charge layering is expected to be a generic 
feature of RTILs at charged interfaces resulting from an 
interplay between steric effects and strong electrostatic 
coupling. While the dependence of the differential ca- 
pacity can be predicted within the mean-field solution 
of the LCG model y, the overscreening and alternate 
charge layering at an electrified interface can not. 

In order to explain overcharging and charge oscilla- 
tions, Bazant et ai. [10] proposed a phenomenological 
theory based on a Landau- Ginzburg-like functional con- 
taining the standard LCG free energy [9| but with an 
additional higher order potential-gradient term, similar 
to what is found in Cahn-Hilliard models, and then solv- 
ing it on a mean- field level. However, this higher order 
gradient term can also be seen [11] as stemming from 
a decomposition of the Coulomb interaction into a long- 
distance mean-field-like component and a non-mean-field 
strong coupling component [12]. This alternative inter- 



pretation motivates a more detailed non-mean-field anal- 
ysis of the original Kornyshev LCG model in order to go 
beyond the limitations of the mean-field approximation. 
We thus propose an exact analysis, albeit in one dimen- 
sion, that demonstrates the full physical phenomenology 
of the LCG model, and in particular shows that over- 
screening emerges naturally from this model without the 
introduction of any new physical interactions. 

We solve exactly the statistical mechanics of a one di- 
mensional LCG in a standard condenser configuration, 
with a method that generalizes other treatments of the 
point-like, or continuous, one dimensional Coulomb gas 
models extensively used in the study of electrolytes [13,]. 
We consider a system where charges q (cations) or —q 
(anions) are located on a line at lattice points with lattice 
spacing a and a total size of M points. The Hamiltonian 
due to Coulomb interactions is 
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where 7 = ^^-^ is the ratio of the electrostatic to ther- 
mal energy. Here Si is a classical spin variable taking the 
value Si = 1 if there is a cation at lattice site z. Si = —1, 
if there is an anion, and 5'^ = if the site is unoccupied. 
We will impose overall electroneutrality on the system, 
stipulating that ^- Si = 0. The grand canonical parti- 
tion function can then be written as 




ms^ 




S = Tr5. 



where /i is the fugacity of both anions and cations, assum- 
ing for simplicity that the ionic liquid is symmetric. Car- 
rying out a Hubbard- Stratonovich transformation, the 
grand canonical partition function can be expressed as 
a path integral over a field (j)j on the lattice, while the 
integral over ip corresponds to the integral over ^0, but 
runs over an interval of length 27r. The partition function 



then assumes the form 
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In the case where there are external charges on the 
boundaries of the system, -\-qQ on the site —1 and —qQ 
on the site M in codenser configuration, we need to add 
i(5(0_i — (/)m) to the action (J4]). The field (j) can be identi- 
fied with the fiuctuating electrostatic potential via the re- 
lation V = —i(p/Pq. One can show that the saddle-point 
equation for the above functional integral, ^I^mf = O5 
reduces to the mean-field equations of Kornyshev [2| and 
Borukhov et al. [9], if one takes the continuous limit 
(a ^ with the maximal charge density q/a kept con- 
stant). 

Writing i/i = (})i and defining p^^'^{y^y') = 
-^= exp ( — ^^~^ ^ J we introduce the symmetric oper- 
ator 

K{y,y')= fp'/^{y,z){l + 2^icos{z))p'/\z,y')dz, (5) 



which allows us to write the grand potential Qq at fixed 
external charges ±Q as 

eM-P^Q) = dx dy e'^^ p^^^K^p^^^ {x,y)e-'^y . 



Introducing the ket- vector I'^q) = p^/^|e 
write 



(6) 
"3!/), we can 



exp(-/?nQ) = (V'Q|if"|^Q). 



(7) 



If instead of fixing the charge at the boundaries one 
keeps fixed the potential difference between them, AF = 
Az///3g, it is easy to show that 



exp(-/31]Ar.) = / c^crexp(-Ai/cr-/3^^) (8) 

where I^az^ is now the grand potential at imposed ex- 
ternal potential difference. It follows that the average 
charge and differential capacity both as a function of the 
imposed potential difference are given by 
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and 



CAu 



d{Q)A. 
dAu 



(9) 



To obtain the charge density we need to know the av- 
erage charge number (Si) on site i. Therefore we just 
need to replace l-\-2ficos{yi) with i 2/isin(?/^) in Eq. (J5]), 
via a new operator 

/CX) 
p'/^{y,z)2fxsin{z)p'/\z,y')dz. (10) 
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FIG. 1: Dimensionless free enthalpy PQm as a function of the 
system size M, for the fugacity /i = 100: solid line: exact, 
Q = 5; dashed line: mean field, Q = 5; dash-dotted line: 
exact, Q = 0; dotted line: mean field, (5 = 0. 



Putting charges Q and — Q on the left and right bound- 
aries, the mean charge density on site i is then 



(5i 



{xI^q\K'LK'^-'-^\xI^q) 

(V'Qliir^lVQ) 



(11) 



In order to compute the thermodynamic quantities de- 
rived here, we need to evaluate the indicated matrix ele- 
ments numerically. The domain of definition of the oper- 
ators K and L depends on the choice of the surface charge 
Q. If we define the Q-dependent class of functions 



keZ-Q 

they are stable under the action of K: 
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The action of the operator K and L, can thus be carried 
out numerically in this Fourier representation. Numeri- 
cal computations are performed by truncating the Fourier 
components at a maximum wave vector /cmax- This ap- 
proximation is valid if 7/c^ax ^ ^ which means that in 
the cases studied here we can take /cmax = 25. 

Numerically one can compute (i) the grand potential 
of the system and the disjoining pressure giving the effec- 
tive interaction between the two bounding layers (ii) the 
differential capacitance as a function of the voltage dif- 
ference between the boundaries and (iii) the dependence 
of the charge density on the lattice position. 

First, we look at the grand potential (^ as a function 
of the size of the system M for given charges iQ on the 
boundaries. The discrete pressure on the boundaries is 
P = Qm — ^M+i and tends to the bulk pressure P5 when 



M is large. The actual force acting on the boundaries of 
the system is given by the disjoining pressure Pd = P—Pb 
that follows from the free enthalpy Qm = ^m + MP^^ 
which is plotted on Fig[Tl at high density and for different 
values of Q. The mean field results are also shown for 
comparison. 

We observe two salient features of the disjoining pres- 
sure: (i) it is attractive when the system size is small and 
goes to zero for M ^2Q (ii) it exhibits oscillations of pos- 
itive/negative pressure when M is even/odd. When the 
size of the system is large, the boundary charges are com- 
pletely screened by the ions, and the plates only feel the 
bulk pressure; when the boundary layers are close to each 
other, not enough space is available for ions to screen the 
boundary charges and they thus attract each other be- 
cause of their opposite charge. Having an even number 
of sites obviously stabilizes the system compared to an 
odd number, however, this effect vanishes as the size in- 
creases. This behavior is not seen within the mean-field 
approximation and is due to the fact that an odd num- 
ber of sites cannot all be filled due to the electroneutral- 
ity condition, and this induces a high energy cost when 
/i is large. When the system becomes large this effect 
is asymptotically negligible, and it does not show up at 
lower densities. The mean field result does not show these 
oscillations, but the general trend agrees with the exact 
result as well as with the continuous model mean-field 
results [14]. Another difference that is not shown here is 
that the bulk pressure is lower in the exact result than 
in the mean-field result. 

The bulk pressure too has some subtle properties. It 
is given by P5 = In(Ao), where Aq is the highest eigen- 
value of K. But the operator K acts on the set of func- 
tions (p!2|) . which depends on the non integer part of Q, 
= Q — [Q\ G [0, 1). Thus we may expect that the bulk 
pressure depends on 0. This is indeed the case, but this 
dependence is very small: AP5/P5 ~ 10~^. This depen- 
dence corresponds to the so-called 6-vacuum introduced 
in Ref. [15| for a standard Coulomb gas. 

We next analyze the differential capacitance as a func- 
tion of the boundary voltage difference for M = 100 sites 
and 7 = 1 at high density (/i = 10) on Fig. [2]and low den- 
sity (/i = 0.1) on Fig. [3l Two phenomena emerge in the 
behavior of capacitance: (i) the capacitance has a dip at 
PZC (ii) the capacitance exhibits oscillations, distinctly 
visible at low densities. The PZC dip appears at low ja 
and confirms the continuum mean-field results of Korny- 
shev [2] . The oscillations stem from the ^-dependence of 
the bulk pressure inducing an extensive ^-dependence of 
the grand potential, which then exhibits minima at 0^ 
corresponding to boundary charges Q G Z + ^*. When 
the imposed voltage varies, the average charge exhibits 
plateaus at these selected charges, and the jumps between 
plateaus induce the peaks in the capacitance. These 
peaks become weaker as the system gets smaller or as 
the temperature increases. 
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FIG. 2: Capacitance as a function of the voltage drop for 
/i = 10 and 7=1: exact result (solid line) and mean field 
result (dashed line). 




FIG. 3: Capacitance as a function of the voltage drop for /j = 
0.1 and 7=1: exact result (solid line) and mean field result 
(dashed line). For smaller 7 the non-monotonicy gradually 
disappears and the exact solution approaches the mean-field 
result of Kornyshev [21 • 



The mean-field capacitance exhibits exactly the same 
trends but without oscillations, which is consistent with 
the absence of plateaus in the average charge. Again our 
mean-field results coincide with those of ref. [2] . 

Next we analyze the charge density profile in the vicin- 
ity of one of the boundaries (left); the same profile with 
opposite charge is obviously found close to the other 
boundary. First we note that for small packing fractions 
(/i = 1) the charge density shows a monotonic variation 
as a function of the separation from the boundary. Fig. 
m and agrees with the mean-field result [9|. 

As we increase the packing fraction (/i > 1, Fig. [5]) 
the overscreening effect starts to dominate the behavior 
of the charge density and clear charge layering emerges: 
a counterion layer is followed by a colon layer with their 




FIG. 4: Mean charge density close to the left electrode (lo- 
cated at X = —1) as a function of the position for Q — 0.5, 
7 = 1 and jjl — 1: exact result (solid line) and mean field 
result (dashed line). 
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FIG. 5: Exact result for the mean charge density close to the 
left electrode (located at x = —1) as a function of the position 
for Q — 0.5, for ^ — 1 and different values of the fugacity. 



thicknesses equal to the lattice unit. This corresponds 
closely to the experimentally observed situation [5] and 
approximately to the results of the model of Bazant et 
al. [10| if the packing fraction is not too large. 

The charge distribution exhibits the following phe- 
nomenology: (i) the amplitude of the oscillations varies 
continuously with the boundary charge Q, being maximal 
for half integer values and zero for integer charges, due to 
the fact that in ID one can have perfect screening only 
for integer values of Q (ii) the overcharging oscillations 
are damped with a characteristic length ^, with ^ ^ ja. 
The characteristic lengthscale of the overscreening oscil- 
lations can be obtained from the two highest eigenvalues 
of the transfer matrix K, Aq and Ai, and the definition 

of the correlation length ^ = ( In ^ j . Analytic com- 



putations of these eigenvalues for /cmax = 1 and (5 = 0, 
and for /cmax = 0.5 and Q = 0.5 indeed gives ^ ~ /i. 

Since the mean-field results show no overscreening, it 
must be a specific feature of the exact result for the same 
LOG model. We thus conclude that in order to describe 
the overscreening phenomena in ionic liquids one needs to 
go beyond the mean-field level of approximation. There 
appears to be no need to modify the original LOG Hamil- 
tonian but it is essential to retain its discretized lattice 
form as opposed to its continuum limit. Based on our 
results, valid in ID, we conclude that an approximation 
akin to the strong-coupling limit, as defined for ordinary 
Coulomb fluids [16], could be pertinent. However, de- 
riving such an approximation systematically for a lattice 
Coulomb gas, where there are multiple length scales, is 
an open problem. 
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